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Abstract 


A simple 2D boundary element analysis, suitable for developing cost effective models for 
tapered composite laminates, is presented. Constant stress and displacement elements are 
used. Closed-form fundamental solutions are derived. Numerical results are provided for 
several configurations to illustrate the accuracy of the model. 


Introduction 

One of the major sources of failure in composite structures is delamination initiating from 
stress concentration sites such as ply drop or a matrix crack. In order to design a damage 
tolerant structure, many parameters affecting delamination must be considered. A 
stacking sequence with the most favorable distribution of interlaminar stresses under 
specified loading, needs to be selected out of a large number of candidate configurations. 
A significant component to achieving affordable rotorcraft composite structures is the 
development of simple and accurate analytical tools that provide trend information at the 
preliminary design stages. 

The analysis of composite structures generally requires numerical modeling. The cost 
associated with conventional finite element simulations of a large number of candidate 
configurations results in a need for developing alternative cost effective models. The 
available engineering beam and shell theories do not allow for a reliable prediction of a 
key parameters such as the peel stress due to simplifying assumptions restricting the 
stress or strain state. Existing approximate elasticity closed-form models do not 
accommodate laminates with tapered geometry. One important application of such 
laminates is composite flexbeams in hingeless and bearingless rotors. During flight, the 
rotor hub arm experiences centrifugal loads as well as bending in the flapping-flexure 
regionl In order to accommodate this bending, the stiffness of the flapping-flexure 
region is changed by varying the thickness of the hub arm. This thickness change is 
accomplished by dropping internal plies in that region. A 2D finite element analysis of 
flexbeam laminates^ requires a mesh with several thousand elements and associated 
number of degrees of freedom of order 10 . 

Cost effective models for elastic analysis of tapered composite structures can be 
developed by applying boundary element techniques. The boundary integral equations 
represent a closed-form model that does not require additional differentiation to obtain a 
solution. Therefore, a coarse mesh and low order elements resulting m small systems of 
linear algebraic equations are expected to provide accurate predictions. Boundaty 
element modeling involves discretization of the domain boundary only. This reduces the 
model size compared to an interior discretization at the same level. On the other hand, 
singular kernel functions have to be used to solve the boundary integral equations^ 
However, this singular formulation is consistent with singular stress fields at the crack 
front of dissimilar ply interface. 
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The first boundary element solutions for plane homogeneous orthotropic elasticity 
problems are provided in References 3 and 4. An analysis for laminated composites in a 
state of strain independent of the longitudinal direction was published in Ref. 5 for the 
case of axial extension and in Ref. 6 for the bending case. Only laminates with straight 
edges were considered. A mesh with 44 elements per ply was sufficient to reach a good 
agreement with classical finite difference and finite element predictions for the 
interlaminar stresses. However, neither of these two references provides the order of 
boundary elements used. 

The main purpose of this work is to develop a simple methodology for applying boundary 
element modeling to rotorcraft composite structures. A 2D boundary integral formulation 
for anisotropic elasticity problems is provided. The constant stress and displacement 
elements, the simplest order element to program, are used. The model simplicity and 
accuracy are first illustrated for a cross sectional analysis of an orthotropic flat laminate 
subjected to axial extension. Results for a plane analysis of generic tapered 
configurations are discussed subsequently. 

Analysis 

In this section, details of a 2D boundary element elasticity modeling of composites are 
presented. The basic relationships are derived for a laminate exhibiting a plane 
deformation state. The cross sectional analysis of laminated beams undergoing a 3D 
strain state independent of the longitudinal direction is simitar but more algebraically 
involved and, therefore, only briefly outlined. 

Boundary Integral Equations 

Consider a 2D elastic domain represented by a longitudinal section of a laminate of an 
arbitrary geometry, undergoing a plane deformation state. Assume that the laminate 
consists of sections of homogeneous material such as differently oriented plies or resin 
pockets. A homogeneous anisotropic material sector will be generally referred to as a 
sublaminate -a ply or a group of plies treated as a single unit with effective properties, in 
the sequel. Equivalent plane elastic properties of a sublaminate can be obtained based on 
the classical lamination theory’ and stiffness or compliance tensor transformation 
relations. The appropriate boundary conditions at the laminate boundary and the material 
interfaces are specified. 

For a sublaminate, the engineering strain-displacement relationships for small 
displacement are 


where x and are Cartesian coordinates and subscript commas denote partial derivatives. 
The following compatibility equations can be obtained from eqns (1) 
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The stress components satisfy the following equilibrium equations 
and the constitutive relations are 
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where S and C denote the compliance and stiffness matrices, respectively. 

In order to obtain the governing integral equations, consider the following two 
equilibrium states: the true state 

0.=C6. (5) 

of a sublaminate subjected to the actual boundary conditions and a known fundamental 
state 


<Tj=CEj (6) 

of a sublaminate of the same material and geometry but subjected to a different set of 
boundary conditions. The only restriction is that the fundamental state has to be singular. 
This condition ensures a non-trivial solution. The fundamental state could be define y 
“cutting” the ply from a half-plane or another simple domain subjected to a concentrated 
force with the point of load application being on the sublaminate boundary. Two 
independent fundamental states associated with the force vector components could be 

obtained. 

The following identity can be obtained from eqns (5) and (6) 

^Ks]z^dxdy = ^<Jlz^dxdy ( 2 ) 

where Cl is the sublaminate area. Substitute eqns (1) into eqn (7), integrate by parts and 
use eqns (3) to obtain 



(8) 
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where r is the sublaminate boundary and 



are the traction and displacement vectors. Equation (8) expresses the reciprocity 
theorem. The simplest possible boundary element model, based on this equation, is 
developed in the following. 

Discretize the ply boundary into a number of segments or elements and assume constant 
ti and U] in each element. Select the middle point of an element as the node at which the 
concentrated load vector is applied to construct the fundamental solution. For each node 
/, two linear algebraic equations can be obtained from eqn (8) 



where Nk is the number of elements, and the superscript indicates the first and second 
fundamental solutions. Summation over the repeated index j is assumed. The total 

M 

number of equations is 2^ TV* where M is the number of sublaminates. The Cauchy 

Jt=l 

principal value of the integrals is implied where necessary. 

There are four parameters per element to be determined: two displacements and two 
traction vector components. However, two of them are specified if the elernent belongs 
to the laminate boundary and four continuity conditions have to be satisfied at the 
sublaminate interfaces. Thus, the boundary conditions reduce the number of unknowns 
to the number of equations. 

It is worth noting that if the tractions and displacements are assumed to change within an 
element according to selected shape functions, eqns (10) have to be appropriately 
modified, and discontinuity points such as comers and crack fronts need to be accounted 

for in the analysis. 

For a laminated beam subjected to a 3D state of strain independent of the longitudinal 
direction, the cross sectional boundary integral equations can be obtained as follows. 

Let Cartesian coordinates x and y be in the cross section plane and z be the longitudinal 
axis. Denote the constant axial strain, bending curvatures and twist rate of the laminate 
by £o > ^ • Assume these four parameters to be zero in a fundamental state. 

Using the same procedure as for the plane deformation, the following reciprocity 
relationship similar to eqn (8) can be derived 
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where the traction and displacement vectors now have three components. The right-hand 
side of eqn (1 1) can easily be transformed to a boundary integral. For an orthotropic 
laminate subjected to axial extension only, with x, y and z coordinate axes parallel to the 
principal material directions 1, 2 and 3, respectively, the boundary element equations 
similar to eqns (10) are 



~ ^0 j" ^\2^2i^y 

r 


where c,y are components of the sublaminate stiffness matrix 
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Fundamental Solutions 

Define the stress function identically satisfying the equilibrium equations (3) as 

CT^ = Fyy ^xy ~ ~^.xy 


(14) 


Substitute the constitutive relations (4), expressed in terms of the stress function (14), 
into the compatibility equation (2) to obtain the following differential equation 

5, - 25„F^, + (25,, + sJF^. - 2s,,F^. + =0 (15) 


The roots of the characteristic equation 
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5 ,,^^ - 25 ,,^^ + ( 25 ,, + 5 „ = 0 


(16) 


are complex® 


=/l,+//y, ^, = ^+i>2 ^3=H ^A=^2 (1”^) 

where the roots (^3 and <^4 are complex conjugate to and . Introduce new variables 

x,=x + A,y y^^M.y * = (18) 

The differential equation (15) can be simplified to 




dxl dyl 


\F = 0 


(19) 


If the characteristic roots and ^2 *^ot repeated, which is the case for anisotropic 
materials, F is a combination of two functions, each being a solution of Laplace’s 
equation, that is 


F = F, + F,, 


d^F d^F, 
dxf "■ dyf 


= 0 


i = \,2 


( 20 ) 


The fundamental solution of Laplace’s equation for a fixed point is a logarithmic function 
of the distance from that point. 

A harmonic stress function Fk is the real part of an analytic function of the complex 
variable z* = Xj, + 1 >* . The derivative of this function 





( 21 ) 


is also an analytic function of z*. The general solution for any plane anisotropic elasticity 
problem can be written in the following form 

<D = 0,(x„y,) + <I>,(x2,y2) (22) 


where 

0,ix„y,) = Ptix„y,) + iQt{x„y,) k = \,2 
are analytic functions. Functions F* and satisfy the Cauchy-Riemann conditions 
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( 24 ) 


dy, dy, dx, 

A fundamental solution for can be written as 

O* = Q In z* = Q In = C, In p, + iC,<p, =P,+ iQ, k = \,l (25) 

where 

Pk =-\l4+yl <p^=arctm— (26) 

The methodology for obtaining the constants in eqns (25) is illustrated for the case of an 
orthotropic material. In addition to the vanishing coefficients and ^25 in the 
compliance matrix (4), the characteristic roots (17) are purely imaginary for a practical 
material system. Equations (18) become 

x^=x y* = Pky (28) 

The following expressions for the displacements and stresses can be derived from eqns 
(1), (4), (14), (18), (20)-(24), and (28) 


u = p^P^ + P2P2 
V = -q\Q\ - qiQi 

^ XX ~ ^\,x ~ ^ ~ ^’2 

(Xy=—Qi.y+—Q2,y 




Pi 


^xy=-P\.y-Pl.y 


where 


(29) 


Pi=-SnP-+Sn ' = 1.2 


(30) 


The first condition to obtain constants C\ and C 2 in eqns (25) is a single-valued vertical 
displacement 

v{(p) = v{(p + 2n) (31) 
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The second condition normalizes the resultant traction force in the horizontal direction in 
the vicinity of the origin of the coordinate system to a half-unit. The resulting 
expressions are 

C=-l C,=-- ^ (32) 

‘ 2 ;t //|92 MiQ2~M2^\ 

Another fundamental solution for is 

o; = -Cli In z, = C>* + iCl In -?- = P** + iQl k = 1,2 (33) 

Pk 


Constants C* and Cj are obtained based on the conditions of a single-valued horizontal 
displacement 

u{(p) = u(<p + 27t) (34) 


and a half-unit resultant traction force in the vertical direction in the vicinity of the origin. 
The resulting expressions are 



1 P2 
2 n p, - P2 



I Pi 
2 n p, - P2 


(35) 


The expressions for the fundamental solutions are summarized in the following 

W* = P\Pk\ 

~^2Qk2 

<^xxk - ~P\Pk\.x ~ P\Pk2,x ^ ~ 

<^yyk-—Qk\.y-^ Qk2.y 

y-2 

^xyk - ~Pk\.y ~Pk2.y 


(36) 


where 
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The same fundamental solutions can be used for the cross sectional analysis of an 
orthotropic laminate subjected to uniform axial extension and bending. The compliance 
coefficients have to be replaced by the following parameters 


Ai ~ ^2 -^n 

S VI 


■*l3-^23 O _ „ ■^23 
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(38) 
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where are coefficients of the compliance matrix for an orthotropic sublaminate. The 

fundamental solutions for the case of general laminates subjected to a 3D strain state 
independent of the longitudinal direction are of a similar form as for the plane case. The 
boundary value problem is governed by two stress functions . The differential equations 
are transformed into three Laplace’s equations and three fundamental solutions have to be 

obtained. 


Results 

In the following, the boundary element model developed in the previous section is 
applied to simple laminate configurations. Numerical results for the interlaminar stresses 
are provided. 

In order to compare the model predictions with published finite element results, consider 
a flat [90/0]s graphite/epoxy laminate subjected to axial extension. The laminate cross 
section is shown in Fig. 1 and the material properties are provided in Table 1 (Ref 9). A 
quarter of the laminate was modeled due to symmetry. The tractions are assumed zero at 
the free edges, the normal displacement and shear stress are set to zero at the symmetry 
axes, and the continuity conditions for both displacements and tractions are satisfied at 
the ply interface. The boundary is discretized as shown in Fig. 2. In order to implicitly 
account for singularities, the nodes are moved at a small distance outside the elements. 
The displacements and tractions are assumed constant within an element as mentioned 
above. As a consequence, all integrals in eqns (12) can easily be calculated analytically. 

The interlaminar stress predictions are compared with results from a hybrid finite element 
cross sectional model^ in Figures 3 and 4. A boundary element mesh that has only 4 
elements per ply at the ply interface, with a total of 24 elements and 48 degrees of 
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freedom, provides accurate predictions at the nodal points. This mesh is labeled coarse 
mesh in Figures 3 and 4, while the fine mesh corresponds to 448 elements with 896 
degrees of freedom. 

Next, a generic tapered glass/epoxy plane configuration shown in Fig. 5 is considered. 
The equivalent plane orthotropic material properties for the ply groups are provided in 
Table 2 (Ref 2). The thick end is clamped and a uniform axial displacement is applied at 
a two tapered length distance from this end. One half of the laminate is modeled due to 
symmetry. The displacements are zero at the clamped end, the tractions are zero at the 
upper surface, and the vertical displacement and the shear stress are zero at the midplane. 

The interlaminar traction predictions at the 0/45 interface are compared with results from 
ABAQUS finite element analysis in Figures 6 and 7. The ABAQUS models are 
generated using constant plane stress four-node reduced integration continuum elements 
(CPS4R). The finite element mesh sizes are 480 elements with 1098 variables that 
include the degrees of freedom and the Lagrange multipliers, and 7712 elements with 
15972 variables. The boundary element model sizes are 72 elements with 144 degrees of 
freedom, and 448 elements with 896 degrees of freedom. The stress discontinuity 
corresponds to the intersection between the uniform and the tapered sections of the 
laminate. A boundary element mesh with 16 elements per ply at the ply interface, 
resulting in a total of 72 elements with 144 degrees of freedom, is sufficient for accurate 
predictions. The model size can be considerably reduced if a graded mesh is used. 

Conclusion 


A 2D boundary element analysis for composite laminates was presented in this work. 
The constant stress and displacement models, corresponding to the simplest element, with 
a small number of degrees of freedom resulted in accurate predictions of the interlaminar 
stresses for the configurations considered. The results presented in this work suggest that 
the boundary element method could become a basis for developing simple and accurate 
models for the elasticity and fracture mechanics analysis of composite structures of 
arbitrary geometry. The accuracy of the boundary element predictions at the contact 
surfaces could result in the development of local models for calculating the strain energy 
release rate components, which will be incorporated in the finite element modeling tools. 
To this end, a basis for simple and accurate plane and cross sectional boundary element 
modeling of composites has been established. The challenge comes with the 
development of closed form fundamental solutions for a 3D analysis of anisotropic 
structures. Generalization of the Stroh formalism could result in such solutions and is a 
subject of future work. 


References 

1. Makeev, A. and Armanios, E. A., “A Simple Elasticity Solution for Predicting 
Interlaminar Stresses in Laminated Composites,” Journal of the American Helicopter 
Society, Vol. 44, No. 2, April 1999, pp. 94-100. 


11 



2. Murri, G. B., O’Brien, T. K., and Rousseau, C. Q., “ Fatigue Life Methodology for 
Tapered Composite Flexbeam Laminates,” Journal of the American Helicopter 
Society, Vol. 43, No. 2, April 1998, pp. 146-155. 

3. Benjumea, R. and Sikarskie, D. L., ”On the Solution of Plane, Orthotropic Elasticity 
Problems by an Integral Method,” Journal of Applied Mechanics, Vol. 39, No. 3, 
September 1972, pp. 801-808. 

4. Snyder, M. D. and Cruse T. A., ’’Boundary-Integral Equation Analysis of Cracked 
Anisotropic Plates,” International Journal of Fracture, Vol. 11, No. 2, April 1975, 
pp. 315-328. 

5. Davi, G., “Stress Fields in General Composite laminates,” AIAA Journal, Vol. 34, No. 
12, December 1996, pp. 2604-2608. 

6. Davi, G. and Milazzo, A., “Bending Stress fields in Composite laminate Beams by a 
Boundary Integral Formulation,” Computers and Structures, Vol. 71, 1999, pp. 267- 
276. 

7. Vinson, J. R. and Sierakowski, R. L., The Behavior of Structures Composed of 
Composite Materials, Martinus Nijhoff Publishers, 1986. 

8. Lekhnitski, S. G., Theory of Elasticity of an Anisotropic Body, Mir Publishers, 1981. 

9. Spilker, R. L. and Chou, S. C., “Edge Effects in Symmetric Composite laminates: 
Importance of Satisfying the Traction-Free Edge Condition,” Journal Composite 
Materials, Vol. 14, No. 1, January 1980, pp. 2-20. 


12 



Table 1. Properties of graphite/epoxy material system’ 


£33 = 20.0 Msi (137.9 GPa) 

£11 =^22 = 2.1 Msi (14.5 GPa) 

G ,3 = G 23 = Gi 2 = 0.85 Msi (5.9 GPa) 


Vj\= Vj2 = Vi2 = 0-21 


Table 2. Orthotropic material properties^ 


Material 

Layup 

£ 11 , GPa 

£ 22 , GPa 

G\ 2 'i GPa 

V12 

S2/E7T1 tape 


47.6 

12.6 

4.81 

0.28 

E-glass/E7Tl-2 fabric 

IBSlHi 

25.3 

24.1 1 

4.56 ^ 

0.153 
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Fig. 1. Flat laminate configuration 
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Fig. 6. Interlaminar surface traction predictions for tapered laminate 



Fig. 7. Interlaminar surface traction predictions for tapered laminate 
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